Effect of strain on surface diffusion in semiconductor heteroepitaxy 
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We present a first-principles analysis of the strain renormalization of the cation diffusivity on the 
GaAs(OOl) surface. For the example of In/GaAs(001)-c(4 x 4) it is shown that the binding of In is 
increased when the substrate lattice is expanded. The diffusion barrier AE(e) has a non-monotonic 
strain dependence with a maximum at compressive strain values (e < 0), while being a decreasing 
function for any tensile strain (e > 0) studied. We discuss the consequences of spatial variations of 
both the binding energy and the diffusion barrier of an adatom caused by the strain field around a 
heteroepitaxial island. For a simplified geometry, we evaluate the speed of growth of two coherently 
strained islands on the GaAs(OOl) surface and identify a growth regime where island sizes tend to 
equalize during growth due to the strain dependence of surface diffusion. 
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I. INTRODUCTION 

In recent years, the heteroepitaxial growth of lattice- 
mismatched semiconductor systems has attracted sub- 
stantial interest. For a number of systems (e.g. Ge/Si, 
InAs/GaAs), it could be shown that under appropri- 
ate experimental conditions the deposited material forms 
elastically strained, dislocation-free (coherent) three- 
dimensional (3D) islands.Em In addition, these islands of- 
ten show a rather narrow size distribution, in particular 
in the higher layers of a stacked 3D array of islands ob- 
tained from repeated deposition of heteroepitaxial films 
separated by spacer layers. This feature is essential for 
the usefulness of these nanostructures as quantum dots 
(QDs) and for their envisaged application in future opto- 
electronic devicesJ3 Considerable theoretical efforts have 
been made in order to rationalize the observed regulari- 
ties in island sizes and ordering. Some approaches have 
attempted to describe the islands as equilibrium struc- 
turesDQ As an alternative explanation, the role of ki- 
netics for the growth of heteroepitaxial islands has been 
emphasized.LrEI It is possible that intrinsic features of the 
kinetics of the growth process give rise to regular struc- 
tures. For instance, self-limiting effects in strained is- 
land growth could result in a preferred-island size, either 
due to a limitation in material supply|13 or due to nucle-. 
ation barriers in the growth of the islands' side facets.til 
This perspective motivated an intense theoretical effort 
towards better understanding the underlying microscopic 
processes in molecular beam epitaxy (MBE), e.g. depo- 
sition, diffusion and nucleation. First-principles calcu- 
lations have already been applied to study different as- 
pects of cation diffusion on unstrained compound semi- 
conductor surfaces. Ej'LJ Up to now, however, the impact 
of strain on the diffusion process still remains elusive, al- 
though attempts-.to make its effect clearer date back to 
the last decade .Ej 

We illustrate the importance of strain for typical het- 
eroepitaxial systems (e.g. Ge/Si, InAs/GaAs) for two 
situations: First, during growth of free-standing het- 



eroepitaxial islands, the islands themselves are under 
compressive strain, whereas the substrate beneath the 
island is expanded. As a consequence of this expansion, 
the substrate surface around an island is under com- 
pressive strain (see, e.g. Ref. fTo]) . Supply of further 
material to the growing island is governed by diffusion 
through this ring-shaped area of compressive strain. A 
suppressive effect of strain on diffusion could slow down 
the growth of larger islands. Secondly, in the growth of 
three-dimensional (3D) stacked arrays of QDs, the buried 
islands act as stressors causing tensile .strain on the cap- 
ping layer in the regions above themJISEZI Again, this may 
affect the growth kinetics of the next layer of islands to 
be formed on the capping layer. Thus, in different stages 
of growth of a nanostructure, different regimes of strain 
may come into focus. 

From the limited published data about strain effects 
on diffusion, it appears surprising that compressive sur- 
face strain could lead to a self-limiting effect on the is- 
land growth. First-principles calculations for diffusion on 
close-packed metal surfaces, in particular Ag/Ag(lll), 
Ref. [l8], have demonstrated that compressive strain in- 
creases the adatom diffusivity by reducing the diffusion 
barrier. Schroeder and WolfEj have extended this finding 
to diffusion on (001) surfaces of simple cubic, fee, and bec 
lattices. Recent molecular dynamics (MD) simulations 
using empirical potentials showed the same trend for Ga, 
In, and As adatompdiffusion on a (2 x l)-reconstructed 
GaAs(OOl) surface. Ej These results also comply with an 
earlier study of Ga kinetics on the strained GaAs(OOl) 
surface ii^ A different strain deaendence of diffusion was 
found, however, for Si adatomEilt3 and dimeia diffusion 
on the Si(001) surface, where tensile strain leads to an 
overall decrease in the diffusion barriers. Yet, the majori- 
ties of the theoretical studies on semiconductor systems 
provide only scarce quantitative information about the 
influence of strain on the diffusion process. 

The aim of this article is two-fold: First, we report 
the results of density-functional theory (DFT) calcula- 
tions for the tracer diffusionEEl of a single In adatom 
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on a GaAs surfaces. In particular, we investigate the 
strain dependence of diffusion in order to clarify the is- 
sues raised above concerning heteroepitaxy of a strained 
system. This problem can be viewed as a 2D analog to 
the effect of pressure on diffusion in bulk materials HJ In 
a second part, we discuss the impact of these findings for 
growth for two typical situations, nucleation on a strained 
capping layer for low In concentrations, and diffusion- 
limited growth of quasi-one-dimensional islands. While 
these topics have been discussed in the4itei'|ature in the 
context of thermo-chemical diffusion,t3n3't2lE3 our focus 
will be on a kinetic description inspired by the results of 
our atomistic calculations. 

For a systematic first-principles investigation of the ef- 
fect of strain on the diffusivity of an In adatem, we de- 
cided to use the GaAs(001)-c(4 x 4) surfaceca as a spe- 
cific example. On top of a complete As surface layer, the 
c(4 x 4) reconstruction has rows of As dimers running in 
the [110] direction, with units of three As dimers inter- 
rupted by a dimer vacancy (see Fig. [I] below) . We have 
chosen this reconstruction because it forms the substrate 
for the initial stages of InAs deposition for temperatures 
T < 500°C (see Ref. |3C]). For very arsenic-rich growth 
conditions, In deposition is expected to lead to direct for- 
mation of 3D InAs islands, as has been shown by previ- 
ous calculations.^ Thus, we can use this system to study 
both the diffusion of the first In atoms on a strained sub- 
strate (e.g. a capping layer with buried islands), as well 
as diffusion in the vicinity of an InAs island on the sur- 
face. We note, however, that the commonly used growth 
conditions to fabricate quantum dots involve formation of 
a InAs wetting layer with reconstructions different from 
the c(4 x 4). Diffusion of In on this wetting layer will be 
addressed in a future publication. 

In Sec. H we outline the underlying computational 
method. The mapping of the potential energy surface 
(PES) for In diffusion on the u nstrai ned Ga As(OO l)- 
c(4x 4) surface is presented in Sec. Ill A. In Sec. IIIB we 
discuss in great details the In adatom interaction with the 
surface As dimers. The effect of strain is then addressed 



in Sec. Ill C . All microscopic results are critically exam- 
ined in Sec. IV in order to assess possible morphological 



consequences for the growth of strained InAs islands on 
GaAs. Finally, a summary and discussion of the results 
is presented in Sec. |v|. 



II. COMPUTATIONAL DETAILS 

DFT calculations have proven to be an efficient tool 
to explore the elementary processes of crystal gmwth 
(see for example RuggeroneJlatsch and SchefflcrE3) . In 
the setting employed here,EJ the substrate is modeled 
by a slab, representing the topmost seven atomic lay- 
ers of the GaAs(001)-c(4 x 4) surface, the bottom layer 
of which was passivated by pseudo-hydrogen atoms. A 
plane-wave basis set with E cut — 10 Ry energy cut- 



off wasp- used in conjunction with ab initio pseudopo- 
tentialsa and the Perdew-Burke-Ernzerhof generalized 
gradient approximation^! to the exchange and correla- 
tion was employed throughout this study. The integra- 
tion over the surface Brillouin zone (SBZ) was performed 
using a Monckhorst-Pack set with two special k-points, 
[(^,0,0), (0, 5,0)], equivalent to 64 k-points in the irre- 
ducible part of the lxl SBZ. Thus, the k-mesh conforms 
withJJie one previously used for the /?2(2 x 4) reconstruc- 
tion.E3 

It is common practice, when addressing adatom dif- 
fusion, to map out the relevant potential energy sur- 
face^) that contains complete information about the dif- 
fusion process. Kley, Ruggerone and SchefflertHl have re- 
cently argued that the adatom-dimer interaction on the 
GaAs(OOl) surface is a crucial factor for the proper de- 
termination of the PES. In order to set the stage for in- 
troducing strain into the problem, we have scanned the 
potential-energy landscape seen by the In adatom, by 
relaxing the latter along the surface normal, placing it 
laterally over a set of equidistant grid points in the (001) 
plane and allowing the topmost six slab layers to freely 
relax. As a reference for calculating the binding energy 
of the adatom we have used the sum of the total energy 
of the (properly strained) bare-surface and the energy of 
a free, spin-polarized In atomHj Geometries were consid- 
ered converged when all residual forces were smaller than 
0.025 eV/A. 

The adatom-dimer interaction poses a multidimen- 
sional problem, since not only the adatom itself, but 
also all degrees of freedom of the surface atoms are in- 
volved in this processes. Therefore, even a full relax- 
ation starting from an adatom above the surface may 
only lead to a local minimum, while other minima may 
exist that can only be reached from different starting con- 
figurations. Test calculations showed that the adatom- 
surface distance and the As- As distance in the dimer are 
most important, and the 2D configurational space defined 
by these coordinates is suitable to image the In-surface 
dimer interaction. Towards this end a special constrained 
relaxation was carried out allowing the In adatom and the 
central As dimer beneath it to be moved as a rigid unit 
(see the inset in Fig. || below). The relative position of 
these three atoms defines a point in a 2D slice through 
the corresponding multidimensional energy hypersurface. 
Performing the constrained relaxation, we succeeded in 
mapping out the 2D PES governing the In-surface inter- 
action in a point-by-point fashion. Further details will 



be given in Sec. [II B 



The strain field at a step or island edge varies slowly 
on the scale of the lattice constant. Therefore we can 
treat the inhomogeneously strained surface by perform- 
ing DFT calculations with a locally adjusted lattice pa- 
rameter. Accordingly, to investigate the influence of 
strain on surface diffusion, the lateral lattice constant 
a was uniformly changed in the range of ±8% around 
its value ao calculated for the unstrained material, thus 
defining the isotropic surface strain tensor e a fj = ed a /3, 
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FIG. 1. (a) Potential-energy landscape for an In adatom 
on the GaAs(001)-c(4 x 4) surface. The adatom is relaxed 
from 2 A above the surface. Contour-line spacing is 0.1 eV; 
atomic positions in the clean surface unit cell are indicated for 
atoms in the upper four layers (As: empty circles; Ga: filled 
circles), side views shown in panel (b), where ao is the bulk 
GaAs lattice constant. The dashed squares show the surface 
unit cell, (c) Sketch of the 2D network of sites used in the 
random walk formalism (four unit cells are depicted). 



with e = a/ao — 1, S a /3 being the Kronecker delta, relax- 
ing again the system as already explained and recording 
the In binding energy for the relevant sites on the PES. 
It is important to note that any change of the super-cell 
volume V is accompanied by a change in the quality of 
the plane-wave basis set. To account for this effect we 
have corrected the calculated total energy of the super- 
cell E to t (E C ut , VLaccording to the scaling hypothesis by 
Rignanese et al.E3 



III. INDIUM DIFFUSION ON GaAs(001)-c(4 x 4) 

A. Potential energy surface 

The mapping procedure resulted in the PES shown in 
Fig. [I] and the binding energies of the In adatom at the 
adsorption sites (Aj) and saddle points (Tk) are given 
in Table | The In adatom, like Ga/GaAs(001)-c(4 x 4) 
(cf. Ref. 13] and Table |J), preferentially adsorbs at the 
four-fold coordinated hollow site Ai (the missing dimcr 



TABLE I. Binding energy Ei (eV) of an In adatom at the 
sites on the GaAs(001)-c(4 x 4) surface, denoted in Fig. 
For comparison, Eb of a Ga adatom at the same sites are 
read off the corresponding potential-energy map calculated by 









Site 








Ax 


A 2 


Ti T 2 


T 3 


T 4 


In 


-2.21 


-1.54 


-1.56 -1.44 


-1.27 


-1.17 


Ga 


-3.04 


-2.20 


-2.54 -2.10 


-2.00 


-1.90 



position), where it interacts with the dangling bonds of 
the second layer As atoms. Two other very shallow min- 
ima A2 are located in between the center and the two 
edge dimers. Jumps between the adsorption sites occur 
through four symmetry-inequivalent saddle points Tfc, 
with Ti being lowest in energy, .— . 

Within transition-state theoryc^l the hopping rate be- 
tween sites Ai and A f , crossing saddle point , is given 
by the ratio of the partition functions of the system with 
the adatom at the equilibrium site, Zj, = Z{A.{), and at 
the saddle point, Zk = Z(Tk), 



r/i = 



1 k B T Z k 
2tt h Zi 



(1) 



where h is the Planck's constant. In the conventional case 
of an elastically relaxed surface the Helmholtz free energy 
F = — k B T hi Z is the proper thermodynamic potential. 
Eq. (fil) can be cast into the well known Arrhenius form 



fi 



^ exp{- 



[£(T fe ) - £(A,)]/fc B T}, 



(2) 



where E is the (static) total energy of the system, read 
off the PES in Fig. 0. The prefactor has the form, 

r£> = exp (-AU vih /k B T + AS vih /k B ) , (3) 



fi 



2tt h 



where At/ V ib and ASVib are, respectively, the associated 

fi 



changes in the vibrational energy and entropy. T\/ is 



temperature-independent within the classical harmonic 
approximation. 

In a simplified picture, the In adatom migrates by a 
random walk on a 2D square lattice defined by the Ai 
sites, Fig. |l| (c). However, we account for hops between 
them via both Ti and T2-A2-T2 with rates Tn and Pn, 
respectively. Indeed, once the In adatom has reached the 
A2 site it needs to overcome a barrier E(T2) — E(A-2) 
of only 0.1 eV in order to move towards a neighboring 
Ai site. As E(T 2 ) - E(A 2 ) < 2k B T for typical growth 
temperatures, the adatom is unlikely to equilibrate at the 
shallow well A2 before it escapes. Thus it is justified to 
use a single rate Tn for the whole path T2-A2-T2. 

Effective diffusion coefficients can now be extracted 
by applying-.the continuous-time random walk (CTRW) 
formalism.BEj Thus, it is easily worked out that the 
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diffusion tensor in Cartesian coordinates (x || [110], 
y || [110]) reads 



0.5 



D a /3 



_ (D [1W] 



D 



[Tio] 



= 4ai 



2 /r u + 4r 11 







(4) 



The factor 4 in front of Tn is partly due to the fact 
that there exist two equivalent T2-A2-T2 paths across 
the block of three dimers. Another factor two enters 
because the In adatom travels in [110] direction \pl- 
times longer distance to reach a neighboring Ai site than 
along the path crossing the saddle Ti. Eq. ([|) thus im- 
plies that an isolated In adatom migrates slightly faster 
in [110] direction, across the dimer rows, than along 
the dimer rows in [TlO] direction, with anisotropy ratio 
■D[iio]/-D(Tio] = 1 + 4Tii/rii. The related diffusion bar- 
riers, entering the rates Tu, and Tn are AE = 0.65 eV 

and AE ~ 0.8 eV (cf. Fig. @ (a), and Table |). One gets 

A9 

a rough estimate for the contribution of the Ai < — > Ai 



channel by assuming that r$ and T\\' differ messen- 
tially. Thus, e.g. at T = 450°C, D [lw] /D^ w] would 
exceed unity by about 50%. At sufficiently low tempera- 
tures, however, one should include A2 in the 2D network 
of sites and consider branching of the diffusion pathways 
towards neighboring Ai or A2 sites. Although an ana- 
lytic result for D a p can still be derived in this case within 
the CTRW formalism, the expressions are rather cum- 
bersome and one has to seek for simplifications requiring 
knowledge of all Tfi rates. |— , 

In comparison to the results by LePage et al.t3 (Ta- 
ble @) for Ga/GaAs(001)-c(4 x 4), In appears to diffuse 
on a less corrugated PES. Thus, for example, at the Ai 
site, being the most stable for both In and Ga, the In 
adatom is by w 0.8 eV less bound than Ga. Further- 
more, the c(4x4) PES provides two additional adsorption 
sites for Ga as compared to In: between the edge dimers 
(Ti in Fig. |) as well as in between Ai and the center 
dimer along [110] (cf. Fig. 2 of Ref. ||). At the Ti site, 
which is a stable adsorption site for Ga, the latter is by 
1.0 eV more strongly bound than In. These differences 
can be easily rationalized in terms of the differences in 
the cation- As bond strength in the corresponding binary 
compoundsJGaAs, InAs) and the larger ionic radius R\ n 
of indiumpPJ we note however that part of the differences 
is to be attributed to the use of the local density approx- 
imation in Ref. |l3|. The cohesive energy per cation- As 
pair is lowest for InAs (E l ^ s = 6.20 eV) as compared to 
GaAs (Eg^ s = 6.52 eV) and AlAs {E^ s = 7.56 eV), 
see e.g. Ref. |4l]. The barriers for diffusion of group-Ill 
cations on the GaAs surface follow the trend given by 
the binding energies, as has been also observed in a first- 
principles study of-Ga and AI diffusion on the GaAs(OOl)- 
(32(2 x 4) surface.B 
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FIG. 2. Binding energy of an In adatom interacting with 
the center As dimer as a function of the ^-coordinate of the 
adatom Zi n . The latter is measured from the z-coordinate of 
the center dimer for the bare surface. Arrows indicate the 
order in which the calculations were performed: adsorption 
(o) and desorption (o). 



B. Interaction of Indium with As-As bonds 

Since the c(4 x 4) reconstruction represents a double 
layer of arsenic, of which the top As atoms form As 
dimers (Fig. [j]), the incorporation of In into the cation 
sublattice requires the topmost arsenic layer to be eventu- 
ally replaced by In atoms. One obvious way how this in- 
corporation could occur is by splitting of the As-As bonds 
in a reaction with an In adatom. For an understanding 
of heteroepitaxy, it is therefore important to study such, 
processes. Furthermore, Kley, Ruggerone and SchefHertHl 
have pointed out for Ga/GaAs(001)-/32(2 x 4) that the 
adatom interaction with the surface As dimers has im- 
portant consequences for the cation diffusivity, since the 
broken As-As bond provides a very stable adsorption site 
for Ga. The underlying mechanism has been identified 
to be the replacement of the rather weak surface As-As 
dimer bond by stronger cation- As bonds, cf. Ref. [l2| For 
a valid description of In diffusion by the PES shown in 
Fig. [j], we therefore have to check if reaction of In with 
the As-As bonds can lead to more stable binding sites for 
In than the minima of the PES. 

We first sample E\, as a function of the adatom height 
Zi n above the T4 site, see Fig. ||. We perform a series 
of calculations for various values of Zi n , where in each 
calculation the adatom is kept fixed, while the substrate 
is allowed to freely relax. In subsequent calculations of 
the series, the geometry of the substrate atoms from the 
previous calculation is used as input. We find that the 
outcome of such a series of calculations depends on the 
initial geometry. While a set of data points modeling ad- 
sorption, starting from Z\ n = 3 A above the closed dimer, 
shows an energy minimum at Z\ n ~ 2.7 A, a series of cal- 
culations for desorption, starting from an adatom incor- 
porated in between the As dimer atoms at Z\ n ~ 0.5 A, 
finds a minimum at Zi n ~ 1.3 A. Both minima have 
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FIG. 3. (a) Binding energy of an In adatom interacting 
with the center As dimer as a function of the As- As distance 
d and the In height above the midpoint of the dimer Zi n , as 
indicated in the inset. (b,c) Bonding configuration and the 
valence electron density in the plane containing the adatom 
and the dimer for the two deeper minima of Eb. 

nearly the same depth, Ef, ~ —1.2 eV. The corresponding 
bonding configurations and total valence electron densi- 
ties are shown in Fig. |^ (b) and (c). The discontinuity 
and hysteresis in Eb seen in Fig. || due to the dimer open- 
ing or closing indicates that the information gained from 
the Z\ n coordinate alone is insufficient for building up a 
complete picture of the adatom-dimer interaction. 

Therefore we have subsequently imaged Eb in a 2D 
configurational space including also the As- As distance d. 
In these calculations, both d and the distance between the 
^-coordinate of the adatom and the As- As midpoint are 
held fixed while relaxing the remaining coordinates of the 
system. As a result, it becomes evident that the adatom, 
approaching the dimer, first goes through a minimum of 
Eb with the dimer bond being almost intact (d = 2.56 A), 
Fig. H (b) . Upon further push towards the surface the In 
adatom splits the dimer bond, overcoming a barrier of 
~ 0.35 eV, and stays shortly within the Z\ n channel in 
a second shallower feature of the PES. The formation of 
directed In- As bonds (see Fig. || (c)) gives rise to a third 
minimum of Eb at (Z\ n = 1.3, d = 5.1 A). 




FIG. 4. Different bonding configurations of an In adatom 
interacting with the edge As dimer; atomic positions in the 
topmost two atomic layers are shown (In: shaded circle, As: 
empty circles), (a) In adatom sitting above the closed edge 
dimer, similar to Fig. ^| (b) (configuration corresponding to 
the T3 site). <f) indicates the As-In-As bond angle; (b) In 
adatom splitting the dimer, similar to Fig. ^|(c); (c) In adatom 
splitting the dimer back-bond. 

In a similar way, we analyze the adatom interaction 
with the edge dimers in the c(4 x 4) unit cell, Fig. || Since 
the edge dimer has only one neighbor, the second-layer 
As atoms are expected to relax more efficiently Indeed, 
comparing similar bonding configurations for the In atom 
at the center dimer, Fig. [| (c), and at the edge dimer, 
Fig. U (b), we find Eb = —1.3 eV for the latter, which 
is only slightly lower than Eb(Ts). A third possibility, 
where In attacks the outer back-bond of an edge dimer, 
is shown in Fig. [| (c) . This configuration results in Eb = 
-1.25 eV. 

It is now clear that additional binding sites for In, re- 
lated to broken As-As bonds, are energetically higher 
than adsorption on the PES of Fig. and are there- 
fore not substantially populated in equilibrium. Thus, 
the mechanism operating in the case of Ga/GaAs(001)- 
(32(2 x 4) is absent in the In/GaAs(001)-c(4 x 4) sys- 
tcmj-Jndeed, the more bulky In adatom with an ionic ra- 
diusEII Rj n larger than that of Ga, when inserted into the 
As dimer, introduces substantial elastic distortion of the 
dimer As back-bonds that cannot be energetically com- 
pensated by the gain due to rehybridization. Note that 
even in the case of an open dimer the In adatom resides 
1.3 A above it, Fig. || (c), which implies an As-In-As bond 
angle <j> ~ 125°, while the Ga adatom is incorporated al- 
most collinearly with the two As atomsjij <j) ~ 175°. In 
summary, our analysis justifies the use of a single PES 
(Fig. |l|) for In diffusion in the CTRW formalism. 

C. Effect of strain 

The foregoing discussion allows us to single out the 
main route for the adatom migration: Ai <-— 1 -> Ai. 
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FIG. 5. (a) Binding energy Eb as a function of strain e 
for an In adatom at the Ai and Ti sites; (b) diffusion bar- 
rier AE = Eb(Ti) — Eb(Ai) as a function of e. Full curves 
on both panels represent least-squares polynomial fits to the 
calculated points. 



Hence, the objective in this section is to analyze the in- 
fluence of surface elastic strain on the Tn rate. A non- 
vanishing strain field in the substrate results in a dif- 
ferent equilibrium configuration of the topmost atomic 
layers. Consequently, both the surface phonon spectrum 
and the PES will experience changes affecting in turn 
both the frequency prefactor and activation energy in 
the exponential in Eq. (Q). The net effect of strain is 
thus determined by the interplay between the latter two 
effects. One may expect, however, that the dominant 
contribution comes from variations in the diffusion bar- 
rier AE = Ef,(Ti) — Eh(Ai), for it enters an exponential. 
This motivated us to concentrate mainly on the strain 
renormalization of AE, but our approach allows for the 
influence of to be also incorporated without detailed 
knowledge of its functional dependence on strain. 

For each particular value of e the In adatom is placed 
above the Ai and Ti sites and the same relaxation 
scheme as for mapping the PES, Sec. Ill A, 



is applied 

to obtain the respective binding energy Eh(e). The cal- 
culated values are shown in Fig. ^. Interestingly, we find 
that AE(e) is a monotonically decreasing function for 
any tensile strain (e > 0) employed in the calculations, 
Fig. H (b). To be specific, this behavior has its onset at 
~ 3% compressive strain, where AE reaches a maximum 
of 0.68 eV. Applying larger compressive strain leads to a 
reduction of AE, with the AE(Q) value recovered again 
for e = —0.06. The non-monotonic dependence on strain 
can be rationalized by inspecting the Eh{e) curves, given 
in Fig. |B| (a). While for e < Eh at the adsorption 
site Ai follows a linear law with a slope of —3.8 eV, the 
binding energy at the saddle point Ti contains, although 
small, non-linear terms in strain which do not cancel in 
the evaluation of AE. For an inhomogeneously strained 
sample, the pronounced strain dependence of Eh for both 
the adsorption site and the saddle point will introduce a 
position dependence of AE. This finding complies with 



none of the two limiting-scenarios of changes of AE dis- 
cussed in the literatureJI3 where only either Eh(Ai) or 
Eh(Tk) contributes. We would also like to emphasize 
that the commonly employed linearityoEl for the strain 
dependence of the diffusion barrier S(AE(e)) is not justi- 
fied in the case of In/GaAs(001)-c(4 x 4) as clearly seen 
from Fig. | (b). Thus one needs to go to higher order 
terms in e to adequately describe the observed S(AE(e)) 
behavior. It is also important to point out that strain 
does not change qualitatively the discussion about the 
interaction of the In adatom with As- As bonds, based 
on the PES in Fig. ||. Extensive tests over the entire 
range of strain considered here were carried out for Eh 
of In at the two stable minima, Fig. |^ (b) and (c). We 
found that the binding configuration of Fig. |1 (b) was 
always slightly preferable over the one in Fig. H (c), but 
the strongest binding site for In remains to be Ai. 

Up to now, the reports in the literature about the 
effect of strain on the diffusion barrier are scarce. A 
slight lowering of the diffusion barrier upon tensile strain 
(e > 3.5%) has been reported for Ag self-diffusion- on a 
Ag(lll) slab within the effective-medium thepry.c3 The 
first-principles treatment of the same system£3 has found 
instead a linear increase of the barrier with strain. For 
semiconductors, a lowering of the diffusion barrier upon 
tensile strain has only beep-,reoorted for Si adatom and 
dimer diffusion on Si(001),EirE2l although the underlying 
binding trends inferred from the MD simulationsEU are 
opposite to those shown in Fig. |B| (a). 

Given the strain dependence of the diffusion barrier 
described above, the basic question arises whether diffu- 
sion limitations can be observed in the growth kinetics 
of InAs on GaAs. This would be the case if the adparti- 
clc diffusivity is reduced for relevant material parameters 
and growth conditions. 

As the substrate around an InAs island, e.g. of pyra- 
midal or truncated pyramidal shape, is under compres- 
sive strain, cf. Ref. ^ and Fig. ^| below, an indium 
adatom approaching the island samples the e < branch 
of AE(e). This branch is accurately described by 



S(AE(e)) = SE n 



e < 0. (5) 



Eq. (||) gives the excess diffusion barrier over the one 
for the unstrained surface AE(0), parameterized by the 
maximum excess 5E m&K = 30 meV, and the strain value 
at which it occurs, e max = —3%. On the basis of Eq. (||), 
one can write a rather general expression for the diffusion 
coefficient taking account of the effect of strain, 



D{e) = D a {l + 2e)5T^\e) 



exp 



KT 



(6) 



where Dq = const is the value of D for the unstrained 
surface, ST^ = T^{e)/T^{0) is the reduction or en- 
hancement factor of the attempt frequency T^-°\ and 
2s = Tr e a p is the relative change in the surface area. 
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imation for the lattice vibrations and a fence-constant 
matrix involving only the degrees of frccdomtfl of the In 
adatom. This estimate indicates ~ 70% increase of the 
prefactor for 4% compressive strain. In this case one can- 
not expect more than a few percent maximum reduction 
of D, which would in turn make diffusion limitations for 
the specific example of the GaAs(001)-c(4 x 4) surface 
rather unlikely. 

IV. CONSEQUENCES FOR GROWTH 

It is interesting to discuss the impact of strain on 
the growth kinetics of both 2D and 3D arrays of self- 
assembled strained islands in a more general context. As 
already mentioned in the introduction, for these two sit- 
uations different regimes of strain are realized, and we 
address them separately. 



FIG. 6. The D/Do = 1 isosurface in the 3D parameter 
space (e,<5r (0) ,T). The view point is from the side where 
D/Do < 1. Points beyond the isosurface correspond to en- 
hanced diffusivity D/Do > 1- 



A first estimate of the expected reduction of D within 
the typical temperature range 350-500° C used for InAs 
deposition on the c(4 x 4)-reconstructed GaAs(OOl) sub- 
strateE] can be obtained by inserting Eq. (ph in (H), set- 
ting <yr(°) = 1. The resulting reduction, D[s ma , x )/ D a ~ 
0.6, turns out to be small due to the smallness of SE max . 
As a consequence, changes of the prefactor due to the ef- 
fect of strain on lattice vibrations are equally important 
in determining the strain renormalization of the In diffu- 
sivity on the GaAs(001)-c(4 x 4) surface in the relevant 
temperature regime. 

Although it is possible to obtain <5r(°) from first- 
principles calculations, it is difficult to get an estimate 
that is better than a factor of two with reasonable com- 
putational effort. Hence we consider ST^ as an inde- 
pendent parameter in the following analysis. Thus the 
right-hand side of Eq. (jfy is a function of three parame- 
ters, e, ST^ , and T. The region in parameter space where 
strain-induced growth limitations can be expected is de- 
fined by the requirement 



D{e- ST^,T)/D <1. 



(7) 



Fig. ^ represents the isosurface D/Dq = 1 in the 3D 
parameter space of e, <Jr^, and T. We note that the 
reduction in diffusivity due to positive S(AE), especially 
for lower temperatures, persists even for<5r(°L> 1, i.e. m 
presence of the so-called compensation effectJl3 However, 
for a very strong compensation effect, 5r^°^(e) > ST^ ~ 
2, no decrease in D in the relevant range of e and T values 
can be expected. 

Finally, we have performed DFT calculations to obtain 
an estimate of (— 0.04), using the harmonic approx- 



A. Growth on a capping layer with buried islands 

In analogy to the previous discussion one can also con- 
sider diffusion in the regime of tensile strain, e > 0. This 
situation is pertinent to the heteroepitaxial growth of 3D 
arrays of vertically self-organized QDs. When a 2D sheet 
of QDs is completed and capped by a spacer GaAs layer, 
the GaAs lattice is expanded in the regions above the 
buried InAs QDs. We consider the onset of In deposi- 
tion onto the spacer layer before nucleation of the first 
new islands. A stationary concentration of adatoms on 
the surface builds up by the equilibrium between supply 
from an atomic In beam source and loss due to evapora- 
tion of In. However, the concentration may vary laterally 
along the surface. In the stationary state the local con- 
centration n(r||) is given by local equilibrium, 

n(r\\) = n exp[-J7(r||)/fc B T], 

where J7(rii) is the binding energy of the In adatoms 
at site Ai, and r 1 1 is the coordinate within the surface. 
f/(r||) is a function of local strain, as given by Eb(e(ru)) 
in Fig. H (a). As can be seen from this figure the binding 
strength increases with increasing strain. Thus, the lo- 
cal concentration of adatoms, and hence the nucleation 
probability for a new island, is increased in the region 
above a buried island where the capping layer surface is 
under tensile strain. Our calculations, thus, provide a 
microscopic foundation for the frequently made assump- 
tionc3't3 that it is easier to nucleate an InAs island on 
those regions of the capping surface where the GaAs 
lattice constant is widened up and thus more closely 
matches the InAs lattice constant. 



B. Diffusion limitations in island growth kinetics 

The conditions under which kinetic growth li mitat ions 
can be expected were discussed already in Sec. HI Q* (see 
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FIG. 7. Strain renormalization of D(x), according to 
Eqs. (^|) and (B), versus island size within the ID flat island 
approximation at T = 450° C and absence of compensation 
effect, <5r'°) = 1. Doc refers to the asymptotic value of D 
at infinitely large distance from the island. The tilt angle 
of the island side facets is 9 and h is their height. The is- 
land/substrate system is assumed to be InAs/GaAs. 

Eq. ([?]) and Fig. ||) . It is interesting to illustrate the possi- 
ble consequences of such limitations for the island sizes. 
This is pursued here within the framework of a simpln 
model problem based on the flat island approximation.!!!! 
One might think, for instance, of adatom diffusion to- 
wards the extended edge of a quantum wire. 

As a first step we address the strain renormalization of 
the adparticle diffusivity due to an isolated island. The 
strain field it creates in the underlying substrate surface, 
within the adopted model, has the form 



e(x) = rj In 



Q 2 (x) 



(8) 



The island geometry (height h, width s, and the tilt an- 
gle of the side facets 9) completely determines the coeffi- 
cients in the second-order polynomials Pi(x) and Q2{x), 
whereas the elastic properties of the material system (e.g. 
Poisson ratio, shear modulus of the substrate) enter the 
prcfactor rj. 

To assess quantitatively the role of diffusion limita- 
tions, we insert the numerical values for <5-E max and e max 
obtained in Sec. 1IIC. For three islands of different size 
Fig. fj] shows the spatial dependence of D{x) obtained 
by inserting Eq. (|) into Eq. (§), for 8T^ = 1 and 
T = 450° C. As seen from the short-range behavior of the 
diffusion coefficient, the larger islands can be about 20- 
30% more effective in hindering the adatom migration 
provided that no substantial compensation effect from 
the r(°) prefactor is present. Consequently, as long as the 
islands grow via strain-dominated surface mass transport 
and the last mentioned condition is met, the compressive 
strain field may lead to retarded growth of the larger 
islands. 

Consider now two islands of characteristic size s\ 
and S2, with s\ > S2, separated by a distance L ^> 
si,S2, Fig. H Supply of adparticles to the surface is 




75 



0.04 > 



distance (nm) 



FIG. 8. Strain field e(x) at the substrate surface, according 
to Eq. (Q, for two ID islands of width si and S2 whose edges 
are separated by distance L. The island/substrate system is 
InAs/GaAs. 



maintained by a stationary flux Fq. One may ask then, 
what is the steady-state adparticle density distribution 
n(x) at the surface, and how does it affect the diffusional 
currents of single adatoms toward the islands —j 1 and 
j2? This is a standard problem in kinetics ; Oc3 how- 
ever, we require it to be solved for a spatially varying 
migration potential U(x) due to the presence of strained 
isla nds. Again, we can exploit the results obtained in 
Sec. Ill Cj , identifying U{x) with Eb(e(x)) for the adsorp- 
tion site Ai. From Fig. [| (a) it becomes clear that a 
compressive strain of a few percent significantly weak- 
ens the binding of the In adatom at the Ai site. Thus 
the coherently strained island gives rise to a repulsive 
potential U{x) that amounts to a few tenths of an eV. 
The time-dependent single atom-density n{x, f) satisfies 
a Smoluchowski-type equationoEil that takes explicit ac- 
count of the field of force due to U(x), 



On 
~dt 



d_ 

dx 



dU(x) 



k^T dx 



Fn 



(9) 



In the simplest case, when the island edges act as per- 
fect sinks, i.e. n(0) = n{L) — 0, the stationary solution 
reads as 



n(x) = Fqc fc B T 



xo-xf , 
e b dx , 



D{x') 



(10) 



with < X{) < L being the position between the two is- 
lands where the total adparticle current vanishes j(xo) cx 
-Vn| xo =0, 



x = 



D(x) 



-dx 



D(x) 



dx, 



(11) 



with D(x) = D(x)exp[—U(x)/k s T]. It is now straight- 
forward to obtain the result that relates j\ with ji 



Xq 



Xo 



(12) 
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Without the effect of strain, the adatom density has a 
simpie parabolic profile, 

F 

n o{x) = 777r( 2x o - x)x (13) 

ZLIq 

with its maximum being exactly at the midpoint between 
the two islands, xq = L/2. The strain renormalization of 
diffusion shifts x towards the bigger island, thus reduc- 
ing the particle current reaching this island. This simple 
ID model problem demonstrates that the smaller island 
will grow faster, until xo gets shifted back towards the 
midpoint when the sizes of the two islands have become 
equal. As a consequence, the strain-limited adatom dif- 
fusion will tend to equalize the island sizes by controlling 
the capture areas for the two islands competing for the 
deposited material. We note, however, that this effect 
will be reduced if a compensation from the frequency 
prefactor r' ^ is operative. 



a wetting layer, and thus the In adatoms, in fact, diffuse 
on an InAs wetting layer or on one of mixed (In,Ga)As 
composition. Therefore the In/GaAs(001)-c(4 x 4) sys- 
tem used here is suitable for modeling the arrival of the 
first In atom to the GaAs surface in the initial stages 
of InAs deposition, but the results cannot be taken over 
directly to the later stages of island growth. Future re- 
search will therefore have to consider diffusion on a wet- 
ting layer and possible account of alloyingEa 
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V. SUMMARY AND DISCUSSION 



We have presented the first ab initio analysis of the 
effect of strain on adatom diffusivity in the context of 
the heteroepitaxial growth of InAs QDs on GaAs(OOl). 
In particular, we quantified the strain dependence of the 
diffusion barrier for an indium adatom on the GaAs(OOT)- 
c(4 x 4) surface. The In interaction with the surface As 
dimers was also given due account. A simple ID model 
problem was employed to demonstrate that the strain- 
limited diffusion contains an archetype of self-limiting 
growth of strained islands. We note that the self-limiting 
effect is an intrinsic feature for a system with spatially 
varying diffusion coefficients, as indicated by Eqs. 
and (12), bearing reduction within a certain operative 
range of strain values. For the lattice-mismatched het- 
eroepitaxy of semiconductor nanostructures, the strain 
effect will give rise to a significant repulsive interac- 
tion between a strained island and an adatom diffus- 
ing towards the island. Moreover, the diffusivity will 
be reduced as well, but a very accurate treatment is 
needed to assess its reduction. Our atomistic calculations 
yielded a maximum increase of the In diffusion barrier 
on GaAs(001)-c(4 x 4) of 30 meV for compressive misfit 
strain. Since this value is of the order of k B T, conclu- 
sions about the relevance of this effect would require as 
well an accurate calculation of the prefactor r^ ^ , which 
is beyond our present goal. 

Finally, we would like to comment on other possible 
extensions of the work presented here. The discussion 
presented in this article refers to diffusivity of—a single 
adatom, the so-called tracer diffusion coefficients (some- 
times denoted in the literature as D*). Our microscopic 
results could, however, serve as an input, for example, to 
Monte Carlo simulations in order to shed more light on 
the effect of strain on adatom diffusivity. Experimentally, 
the self-assembled coherent islands are usually grown on 
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